ParamSe=zeros(NZ,1); %Standard error of estimated coefficients of Z variables
MeanSE=zeros(NV,1); %Standard error of Mean
StdSE=zeros(NV,1);  %  " of Std Devs
CorrSE=zeros(NV,NV);  % " of Correlation matrix
FreqSE=zeros(NV,NBins); % " of Freq in each bin
MidSE=zeros(NV,NBins); % " of Midpoint for each bin, (std err should be zero)

paramhold=zeros(NZ,NReps);
llhold=zeros(1,NReps);
mnhold=zeros(NV,NReps);
stdhold=zeros(NV,NReps);
corrhold=zeros(NV,NV,NReps);
freqhold=zeros(NV,NBins,NReps);
midhold=zeros(NV,NBins,NReps);

XMAT_Original=XMAT;
PID_Original=XMAT(:,1);
CSID_Original=XMAT(:,2);
clear i r

llorig = finalLL;
pricetime0 = pricetime;
qualitytime0 = qualitytime;
wlifetime0 = wlifetime;
infratime0 = infratime;
jobstime0 = jobstime;
buffertime0 = buffertime;
swtime0 = swtime;
captime0 = captime;

igpricetime0 = igpricetime;
igqualitytime0 = igqualitytime;
igwlifetime0 = igwlifetime;
iginfratime0 = iginfratime;
igjobstime0 = igjobstime;
igbuffertime0 = igbuffertime;



for xi=1:NReps
    FirstRun=0;
    disp('Bootstrap Sample #');
    disp(xi);
    disp(datetime);
    disp('');
    btsample=randi(NP,[NP,1]);
    XMAT=[];
    
    pricetime = [];
    qualitytime = [];
    wlifetime = [];
    infratime = [];
    jobstime = [];
    buffertime = [];
    swtime= []; 
    captime= [];
    
    igpricetime = [];
    igqualitytime = [];
    igwlifetime = [];
    iginfratime = [];
    igjobstime = [];
    igbuffertime = [];
    
    for r=1:NP;
       thisp=btsample(r,1);
       rowsr = PID_Original==thisp;
       thisx=XMAT_Original(rowsr,:);
       thisx(:,1)=r.*ones(size(thisx,1),1); %creates the ID number sequentially
       XMAT=[XMAT ; thisx];   
       pricetime = [pricetime; pricetime0(rowsr,:)];
       captime = [captime; captime0(rowsr,:)];
       swtime = [swtime; swtime0(rowsr,:)];
       buffertime = [buffertime; buffertime0(rowsr,:)];
       qualitytime = [qualitytime; qualitytime0(rowsr,:)];
       wlifetime = [wlifetime; wlifetime0(rowsr,:)];
       infratime = [infratime; infratime0(rowsr,:)];
       jobstime = [jobstime; jobstime0(rowsr,:)];
       
       igpricetime = [igpricetime; igpricetime0(rowsr,:)];
       igbuffertime = [igbuffertime; igbuffertime0(rowsr,:)];
       igqualitytime = [igqualitytime; igqualitytime0(rowsr,:)];
       igwlifetime = [igwlifetime; igwlifetime0(rowsr,:)];
       iginfratime = [iginfratime; iginfratime0(rowsr,:)];
       igjobstime = [igjobstime; igjobstime0(rowsr,:)];
    end
    XMAT(:,2)=CSID_Original;
    if discountfun=="QuasiHyper"
        CreateDrawsDisWtpProbsQHpref;
        CreateZQH;
    else
        CreateDrawsDisWtpProbs;
        CreateZ;
    end
    if YesGPU==1       
        PROBS=gpuArray(PROBS);
        BETAS=gpuArray(BETAS);
        Z=gpuArray(Z);
    end
    param=paramhat_orig; %StartB; %Trying to start from estimate of real dataset
    options=optimset('LargeScale','off','Display','final','GradObj','on',...
    'MaxFunEvals',10000,'MaxIter',MAXITERS,'TolX',PARAMTOL,'TolFun',LLTOL,'DerivativeCheck','off');
    [paramhat,fval,exitflag]=fminunc(@flexll,param,options);
    if YesGPU==1;
        PROBS=gather(PROBS);
        BETAS=gather(BETAS);
        Z=gather(Z);
    end;
    paramhold(:,xi)=paramhat;
    llhold(1,xi) = finalLL;
    [mnhold(:,xi),stdhold(:,xi),cc,freqhold(:,:,xi),midhold(:,:,xi)]=stats(paramhat,NBins);
    if discountfun=="QuasiHyper"
        corrhold(:,:,xi)=corrcov(cc);%pref space
    else
        corrhold(:,:,xi)=corrcov(cc(3:end,3:end));
    end
    %clear f;
    save(filename);
end
finalLL = llorig;

igqualitytime = igqualitytime0;
igwlifetime = igwlifetime0;
iginfratime = iginfratime0;
igjobstime = igjobstime0;
igbuffertime = igbuffertime0;
igpricetime = igpricetime0;

qualitytime = qualitytime0;
wlifetime = wlifetime0;
infratime = infratime0;
jobstime = jobstime0;
buffertime = buffertime0;
swtime = swtime0;
captime = captime0;
pricetime = pricetime0;

XMAT=XMAT_Original;
clear llorig XMAT_Original qualitytime0 wlifetime0 infratime0 jobstime0 buffertime0 swtime0 captime0;
clear igqualitytime0 igwlifetime0 iginfratime0 igjobstime0 igbuffertime0 igswtime0 igcaptime0;

ParamSE=std(paramhold,0,2);
MeanSE=std(mnhold,0,2);
StdSE=std(stdhold,0,2);
CorrSE=std(corrhold,0,3); 
FreqSE=std(freqhold,0,3); 
MidSE=std(midhold,0,3);

disp('Final Log-likelihood');
disp(' ');
disp(finalLL);
disp(' ');
disp('Utility Parameters');
disp(' ');
disp('                    Mean                   Std Dev');
disp('              ------------------   -----------------------');
disp('                 Est     SE            Est         SE');

for r=1:NV
        fprintf('%-10s %10.4f %10.4f %10.4f %10.4f\n', NAMES{r,1}, [MeanEst(r,1) MeanSE(r,1) StdEst(r,1) StdSE(r,1)]);
end
if discountfun=="QuasiHyper"
    disp(' ');
    disp('Correlations');
    disp(corrcov(CovMatEst));
    disp('T-stats on Correlations');
    disp('Diagonal has Inf because diagonals of correlation matrix are always 1.');
    disp(corrcov(CovMatEst)./CorrSE);
    disp(' ');
else
    disp(' ');
    disp('Correlations');
    disp(corrcov(CovMatEst));
    disp('T-stats on Correlations');
    disp('Diagonal has Inf because diagonals of correlation matrix are always 1.');
    disp(corrcov(CovMatEst)./CorrSE);
    disp(' ');
end
Lm=zeros(NV,1);
Um=zeros(NV,1);
Ls=zeros(NV,1);
Us=zeros(NV,1);
Lp=zeros(NZ,1);
Up=zeros(NZ,1);
for(ii=1:(NZ))
    Lp(ii,1)=quantile(paramhold(ii,:),0.025);
    Up(ii,1)=quantile(paramhold(ii,:),0.975);
end
for(ii=1:NV)
    Lm(ii,1)=quantile(mnhold(ii,:), 0.025);
    Um(ii,1)=quantile(mnhold(ii,:), 0.975);
    Ls(ii,1)=quantile(stdhold(ii,:), 0.025);
    Us(ii,1)=quantile(stdhold(ii,:), 0.975);
end

BootCIparam = horzcat(paramhat_orig, ParamSE, Lp,Up);
tab = dataset({BootCIparam 'Estimate', 'SE', 'Lower', 'Upper'},...
    'obsnames', {});
tab = dataset2table(tab);
nm = horzcat('ParamResults.',modelname, '.csv');
writetable(tab, nm);

%disp('First Row of ParamResults.csv contains Price coefficent if in preference space')

BootCImean = horzcat(MeanEst, MeanSE, Lm,Um, StdEst, StdSE, Ls, Us);

Results = dataset({BootCImean 'MeanEst', 'MeanSE', 'LowerMean', 'UpperMean', 'SDEst', 'SDSE', 'LowerSE', 'UpperSE'},...
    'obsnames', NAMES);

Results = dataset2table(Results);
nm = horzcat('DistrResults.',modelname, '.csv');
writetable(Results, nm,'WriteRowNames',true);

disp(Results);
save(filename)

